library(lmtest)
library(sandwich)
library(stargazer)
library(tidyverse)


#ROBUSTNESS CHECK####
rm(list = ls())
load("figure4.RData")

ds <- list(final_90, final_105, final_120, final_135, final_150, final_165, final_180,
           final_195, final_210, final_225, final_240, final_255, final_270, final_285,
           final_300, final_315, final_330, final_345, final_360)

rdd_coef <- rdd_se <- NULL
rd_times <- seq(90, 360, 15)

for(i in 1:length(ds)){
  model <- lm(afd_second ~ as.factor(condition), as.data.frame(ds[i]))
  coefs <- coeftest(model, vcov = vcovHC(model, type="HC1"))
  rdd_coef <- c(rdd_coef, coefs["as.factor(condition)treatment", 1])
  rdd_se <- c(rdd_se, coefs["as.factor(condition)treatment", 2])
  rm(model, coefs)
}

to_plot <- data.frame(rd_times, rdd_coef, rdd_se)
to_plot$upper <- to_plot$rdd_coef + 1.96*to_plot$rdd_se
to_plot$lower <- to_plot$rdd_coef - 1.96*to_plot$rdd_se


ggplot(to_plot, aes(rd_times, rdd_coef)) +
    geom_point() + 
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "darkblue") + 
    labs(x = "Days since/to elections", y = expression(hat(beta)~" and 95% CI")) +
    theme_bw() +
    theme(axis.text.y =  element_text(size = 12), axis.text.x = element_text(size = 12), 
          axis.title.y =  element_text(size = 14), axis.title.x = element_text(size = 14))+
    scale_x_continuous(breaks = seq(60, 360, 30))
